// Optional MACC graph
// Author Levi Marks levi.marks.1>at>gmail.com

// Value of conserved gas for area plot
set obs `=_N+1'
replace marginal_value = marginal_value[_n-1] if missing(marginal_value)
replace reduction_percent = 100 if missing(reduction_percent)

// Indicating out-of-sample predictions: Only 79% that is methane is taxed, so the highest observe price (5.75) is equivalent to a methane tax of 5.75/.79 = 7.278
set obs `=_N+1'
replace reduction_percent = 97.75 if _n == _N
gen reduction_percent_end = 97.75 if _n == _N
replace mc_mcf = 7.278 if _n == _N
gen mc_mcf_end = 29.75 if _n == _N

set obs `=_N+2' // This makes a little dash
replace mc_mcf = 7.278 if _n >= _N - 1
replace reduction_percent = 97 if _n == _N - 1
replace reduction_percent = 98.5 if _n == _N

twoway area marginal_value reduction_percent, color(gs14) lwidth(vvthin) lcolor(gs14) || ///
	line mc_mcf reduction_percent if _n < _N - 3, lcolor(black) lwidth(vthin) || /// Thin line then thick line to double check axes are scaled correctly
	line mc_tco2e ag_reduction_tco2e, lcolor(gs3) lwidth(medthick) yaxis(2) xaxis(2) || ///
	line mc_mcf reduction_percent if _n >= _N - 1, lcolor(gs12) lwidth(thin) || ///
	pcarrow mc_mcf reduction_percent mc_mcf_end reduction_percent_end if _n == _N - 2, lcolor(gs12) mcolor(gs12) mlwidth(thin) msize(1.75) barbsize(.35) lwidth(thin) ///
	yline(5.68, lcolor(gs8) lpattern(dash)) /// Obtain these by examining mc_mcf in the dataset at i.e. tax_tco2 = 5
	yline(12.67, lcolor(gs8) lpattern(dash)) ///
	yline(25.41, lcolor(gs8) lpattern(dash)) ///
	text(5.68 2 "$5 Carbon Price" , yaxis(1) xaxis(1) placement(ne) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(12.67 2 "$20 Carbon Price", yaxis(1) xaxis(1) placement(ne) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(25.41 2 "Social Cost of Methane", yaxis(1) xaxis(1) placement(ne) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(.5 50 "Value of Captured Gas", yaxis(1) xaxis(1) placement(n) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(`=(29.75+7.278)/2' 97.5 "Out-of-Sample" , yaxis(1) xaxis(1) placement(w) orient(vert) margin(.5 .5 .5 .5) color(gs12) size(*.8)) ///
	graphregion(color(white)) legend(off) xsize(6) ysize(4) ///
	plotregion(margin(0 0 0 0)) ///
	xlabel(0 "0%" 20 "20%" 40 "40%" 60 "60%" 80 "80%" 100 "100%", axis(1) labsize(medsmall) format(%12.0fc) nogrid) ///
	xlabel(15000000(15000000)60000000 , axis(2) labsize(medsmall) format(%12.0fc) nogrid labcolor(gs8) tlcolor(gs8) tlwidth(medthin)) ///
	xscale(range(0 100)) ///
	xscale(range(0 `=total_emissions_tco2e') axis(2) lcolor(gs8) lwidth(medthin)) /// Calculate total CO2e emissions for this one
	xtitle("Total Abatement"  , axis(1) size(medlarge) margin(0 0 -2 2)) ///
	xtitle("Total Abatement (tCO2e)", axis(2) size(medlarge) margin(0 0 2 -2) color(gs8)) ///
	ylabel(0(5)31, axis(1) labsize(medsmall)format(%6.0fc) nogrid) ///
	yscale(range(0 31)) ///
	ylabel(10(10)40, axis(2) labsize(medsmall) format(%6.0fc) nogrid labcolor(gs8) tlcolor(gs8) tlwidth(medthin)) ///
	yscale(range(0 `=31*53.68/34') axis(2) lcolor(gs8) lwidth(medthin)) ///
	ytitle("Marginal Cost ($/Mcf)"  , axis(1) size(medlarge) margin(0 2 0 0)) ///
	ytitle("Marginal Cost ($/tCO2e)", axis(2) size(medlarge) margin(2 0 0 0) color(gs8)) 
